Finite-difference modeling experiments for seismic interferometry

In passive seismic interferometry, new reflection data can be retrieved by crosscorrelating recorded noise data. The quality of the retrieved reflection data is, among others, dependent on the duration and number of passive sources present during the recording time, the source distribution, and the source strength. To investigate these relations we set up several numerical modeling studies. To carry out the modeling in a feasible time, we design a finite-difference algorithm for the simulation of long-duration passive seismic measurements of band-limited noise signatures in the subsurface. Novel features of the algorithm include the modeling of thousands of randomly placed sources during one modeling run. The modeling experiments explore the dependency relation between the retrieved reflections and source-signature length, source positions, number of sources, and source amplitude variations. From these experiments we observed that the positions of the passive sources and the length of the source signals are of direct influence on the quality of the retrieved reflections. Random amplitude variations among source signals do not seem to have a big impact on the retrieved reflections. INTRODUCTION Seismic interferometry (SI) is a relatively new branch of geophysics and deals with the retrieval of new seismic responses between receivers by crosscorrelating responses recorded at these receiver locations. Applications of SI exist for exploration data with controlled sources (Schuster et al., 2004; Bakulin et al., 2004, 2006; Wapenaar, 2006), as well as for passive data due to natural sources (Rickett and Claerbout, 1999; Wapenaar et al., 2002; Draganov et al., 2007). For a cross-discipline overview, we refer the readers to the Geophysics Reprint Series about SI (Wapenaar et al. (editors), 2008), which also contains contributions from authors from other disciplines. The retrieval of surface waves using natural sources has already led to numerous successful studies (e.g. Campillo and Paul [2003]; Shapiro and Campillo [2004]). SI for the retrieval of reflected body waves from passive measurements has been successfully applied only recently (Draganov et al., 2007, 2009). However, in many cases it remains difficult to interpret the retrieved wavefields and verification with modeled or measured results remains useful. To test ideas and new concepts for SI, aiming at the retrieval of reflected body waves, passive measurements are needed. Despite the many advantages that SI could offer, passive measurements are still rare and for the moment, we have to rely on numerical forward-modeling studies to gain experience in the practical use of passive SI. The goal of the modeling studies discussed in this paper, and carried out by the specially designed algorithm, is to get a better understanding of what influences the quality of the retrieved reflections. In the following sections we give insights to the relation between the quality of the retrieved reflections and: • the average duration of the passive sources; • the number of the passive sources captured during the recording time; • the source distribution; • the presence of intrinsic attenuation in the medium; • variations in the passive-source amplitudes; • the effect of receiver topography. This list is by no means complete and represents first steps in studying and quantifying the quality of retrieved reflections. For the simulation of passive measurements, very long recording times are needed (from minutes to hours), many (thousands) of random source positions, random source signatures, and random excitation times have to be included in the modeling algorithm. Without the software, presented as a part of this paper, we would not have been able to perform all the experiments discussed in the following sections. Along with the discussion about the experiments, implementation aspects of the algorithm are explained as well. To guide the reader, Table 1 lists the kind of experiment carried out and the figure number showing the results of the experiment. Manuscript received by the Editor 28 January 2010; revised manuscript received 10 June 2011; published online 28 December 2011. Delft University of Technology, Department of Geotechnology, Delft, Netherlands. E-mail: [email protected]; [email protected]. © 2011 Society of Exploration Geophysicists. All rights reserved. H1 GEOPHYSICS. VOL. 76, NO. 6 (NOVEMBER-DECEMBER 2011); P. H1–H18, 19 FIGS., 1 TABLE. 10.1190/GEO2010-0039.1 Downloaded 26 Jan 2012 to 131.180.61.65. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ We start this paper by briefly reviewing the Green’s function representation for SI with two-way wavefields. The first example in the “Modeling Experiments” section illustrates the basic concepts of SI. We use the result in that example as a reference for comparison with subsequent retrieved SI results. We then continue by presenting the results for different variations of random source positions and source signatures and investigate their connection to the quality of the retrieved reflection data. We will also explain how these experiments can be carried out using the accompanying software. The manual, bundled with the software, should be consulted for more detailed explanations about the parameter settings to use. The last section shows SI results for a more complicated model and includes topography. In the Appendix, the finite-difference algorithm and the generation and implementation of band-limited random-noise signals is explained. SI WITH TWO-WAY WAVEFIELDS In the brief theoretical background given in this section, we show the main relations for SI as derived by Wapenaar and Fokkema (2006) from reciprocity theory. Consider a Green’s function Gðx; xA; tÞ for an inhomogeneous lossless acoustic medium, where x and xA are the Cartesian coordinate vectors for the observation and source points, respectively, and where t denotes time. We define the temporal Fourier transform as Ĝðx; xA;ωÞ 1⁄4 ∫ −∞ expð−jωtÞ Gðx; xA; tÞdt, where j is the imaginary unit and ω the angular frequency. Assuming the unit point source at xA is of the volumeinjection-rate type, the wave equation for acoustic pressure Ĝðx; xA;ωÞ reads ρ∂iðð1∕ρðxÞÞ∂iĜðx; xA;ωÞÞ þ ðω2∕c2ðxÞÞĜðx; xA;ωÞ 1⁄4 −jωρδðx − xAÞ: (1) Here cðxÞ and ρðxÞ are the propagation wavespeed and the mass density of the inhomogeneous medium and ∂i denotes the partial derivative in the xi-direction (Einstein’s summation convention applies to repeated subscripts). In the remainder of this section we will leave out the dependency on frequency in the notation of the Green’s functions Ĝ. The representation of Ĝ, as derived for SI from Rayleigh’s reciprocity theorem (Rayleigh, 1878; Bojarski, 1983; Wapenaar et al., 2004, 2005), reads ĜhðxA; xBÞ 1⁄4 I ∂D −1 jωρðxÞ ðĜ ðxA; xÞ∂iĜðxB; xÞ − ð∂iĜ ðxA; xÞÞĜðxB; xÞÞnidx; (2) with ĜhðxA; xBÞ1⁄4 ̂ ĜðxA; xBÞ þ Ĝ ðxA; xBÞ 1⁄4 2RfĜðxA; xBÞg; (3) where ∂D is an arbitrary closed surface with an outward pointing normal vector n 1⁄4 ðn1; n2; n3Þ and the asterisk denotes complex conjugation. The points xA and xB are both situated inside ∂D; the medium may be inhomogeneous inside as well as outside ∂D. Note that equation 2 is exact and applies to any lossless arbitrary inhomogeneous acoustic medium. The retrieved Green’s function ĜðxA; xBÞ contains, apart from the direct wave between xA and xB, all scattering contributions including primary and multiple reflections from inhomogeneities inside as well as outside ∂D. To make equation 2 more suited for practical applications, several approximations are made: The medium at and outside ∂D is assumed to be homogeneous, with propagation velocity c and mass density ρ; a far-field approximation expresses the dipole sources in terms of scaled monopoles; the medium parameters around ∂D are assumed to be smoothly varying (Wapenaar and Fokkema, 2006). These approximations lead to the following equation 2RfĜðxA; xBÞg ≈ 2 cρ I ∂D Ĝ ðxA; xÞĜðxB; xÞd2x: (4) In general, these approximations involve amplitude errors that can be significant (Ramirez and Weglein, 2009). Furthermore, spurious events may occur due to incomplete cancellation of contributions from different stationary points when the medium outside ∂D is inhomogeneous (Draganov et al., 2004). However, since the approximations do not affect the phase of equation 4, it is considered acceptable for SI. When we assume that the sources are uncorrelated (both in space and in time) we can write the observed wavefields as ûðxAÞ 1⁄4 I ∂D ĜðxA; xÞN̂ðxÞd2x and ûðxBÞ 1⁄4 I ∂D ĜðxB; xÞN̂ðxÞd2x; (5) where the noise spectrum N̂ðx;ωÞ has to fulfill hN̂ðxÞN̂ ðx 0Þi 1⁄4 δðx − x 0ÞŜðωÞ; (6) where h:i denotes a spatial ensemble average and ŜðωÞ the power spectrum of the noise sources. Equation 6 states that N̂ðxÞ and N̂ðx 0Þ are mutually uncorrelated for any x ≠ x 0 at ∂D, and that their power spectrum ŜðωÞ is the same for all x. Using equations 5 and 6 in equation 4 results in 2RfĜðxA; xBÞgŜðωÞ ≈ 2 cρ hûobs ðxAÞûðxBÞi: (7) If part of ∂D is formed by a free surface (∂D0), then in equation 4 and 5 we will need sources only of the part ∂D1 1⁄4 ∂D − ∂D0 that Table 1. The different experiments carried out and the figure number showing the results. Experiment Figure number source duration Figure 4 number of sources Figure 5 and 6 source distribution Figure 7, 8, and 9 attenuation Figure 10 source amplitude Figure 12 and 13 receiver topography Figure 15 H2 Thorbecke and Draganov Downloaded 26 Jan 2012 to 131.180.61.65. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ does not include the free surface (see Figure 1). When the receivers xA and xB are placed on the free surface, equation 7 remains valid and the observed quantities u are then particle velocities (Wapenaar and Fokkema, 2006). We use equations 7 and 5 in the remainder of this paper to retrieve reflection data from passive measurements. MODELING EXPERIMENTS The modeling experiments in this section investigate how retrieved reflections depend on different passive-source configurations. At the same time, it illustrates the capabilities of the finitedifference algorithm bundled with the paper. To get an idea of how well the reflections could be retrieved with passive SI, we would like to compare the results that we obtain with an active SI reference result obtained using a regular distribution of transient sources along a closed contour (shown as a dashed line in Figure 2) around particle-velocity receivers placed at the free surface. A source is placed at every grid point on this dashed line, with a grid spacing of 10 m. The sources are modeled sequentially using a loop to start a new FD modeling for each source position. The source signal is a first derivative of a Gaussian wavelet with a frequency peak at 13 Hz. The receivers are placed on the free surface at z 1⁄4 0 on a 50 m grid covering the whole surface. This experiment follows the theory as close as possible according to interferometry equation 4. Figure 3 shows the contribution of the different parts of the closed contour to the retrieved reflections. Figure 3a shows the contribution of sources placed on the lower horizontal part (z 1⁄4 3600 m) of the contour. This result shows that the higher angles in the reflection response are missing, but all reflections and multiples are already present. The contribution of the two vertical parts of the closed contour in Figure 3b mainly retrieves the higher angles of the reflections. Outside and close to the upper part of the vertical edges of the chosen contour the assumptions in equation 4 are violated, and spurious events can be observed. In Figure 3b these spurious events are clearly visible outside the integration contour, between −5000∶ − 4000 and 4000∶5000 m and exhibit a reverse curvature of the reflection arrivals. The contribution of the complete closed contour is shown in Figure 3c and we can see, as expected, a successful retrieval of all reflection events. For comparison, in Figure 3d a directly forward-modeled result is shown. Influence of the number of sources and source-signature length As a first SI experiment with passive sources, we use again the model shown in Figure 2 but with random source positions below level z 1⁄4 500 m. In Figure 2 the source positions, in total 1000, are shown as black dots. In the algorithm a square region, where the sources can be placed, is defined by the four corners using the parameters xscr1, xsrc2, zsrc1, zsrc2. For the investigation of the sources’ influence on the retrieved result, the source signal duration and start time is varied, while the source strength is the same for all sources. The source signature is a random sequence with a maximum frequency of 30 Hz and constructed according to the procedure explained in the Appendix. The FD program Figure 1. Green’s function (Ĝ) retrieval by crosscorrelation requires sources (★) on a closed surface. When part of the closed surface is a free surface (∂D0), it suffices to have sources on the remaining part (∂D1). The rays represent the full responses, including primary and multiple scattering due to inhomogeneities and reflections from the free surface ∂D0. -500 0 50